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ABSTRACT 


This  report  explores  the  algebraic  and  measure-theoretic 
properties  of  the  Euler  parameter  (quaternion)  representation  of 
the  rotation  group.  A  family  of  probability  densities  on  this 
group  is  examined,  and  schemes  for  numerical  integration  and  esti¬ 
mation  of  satellite  attitude  are  presented  with  numerical  examples. 
An  effort  has  been  made  wherever  possible  to  deal  with  nonlinear 
problems  directly  rather  than  linearizing  them.  Several  questions 
for  further  research  are  raised  regarding  convolutions  of  densities 
on  the  group.  Brownian  motion,  and  the  statistics  of  spin  and  torque 
vectors  considered  as  stochastic  variables. 
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1.  INTRODUCTION 


Euler's  Theorem  states  that  the  most  general  rotational 
(or  angular)  displacement  of  a  rigid  body  can  be  accomplished  by  a 
single  rotation  through  an  angle,  9,  about  a  line.  Letting  the 
direction  of  this  line  be  designated  by  a  unit  vector 


we  have  a  numerical  description  of  this  rotational  displacement 
(henceforth,  just  "displacement").  It  is  not  unique,  however, 
since  both  u  and  -u  define  the  same  line  and  0  has  not  been  re¬ 
stricted  to  an  interval  of  length  2ir.  This  representation  is  not 
particularly  economical,  using  four  variables  and  one  constraint, 

]u|  =1,  to  describe  a  mathematical  object  that  has  three  degrees 
of  freedom. 

There  are  many  choices  for  parameterizing  the  rotation 
group  with  three  variables.  One  of  the  most  common  is  the  z-x-z 
Euler  angles,  three  plane  rotations  specified  in  "running  coordi¬ 
nates;"  another  is  the  roll-pitch-yaw  system,  or  x-y-z  Euler  angles, 
used  in  models  of  ship  motion.  These,  too,  have  their  drawbacks 
stemming  from  the  differences  between  the  topology  of  this  group 
3 

and  that  of  R  ,  the  real  3-dimensional  space  of  the  three  angles. 
Such  representations  cannot  be  uniformly  continuous  —  "small" 
rotations  are  not  always  described  by  small  parameter  values.  This 
presents  numerical  difficulties  when  values  of  these  parameters  are 
being  estimated  from  observations  of  the  orientation  of  a  rigid 
body.  Furthermore,  the  propagation  equations  for  rigid-body  motion 
are  not  as  simple  as  one  would  like  and  suffer  from  the  same  topo¬ 
logical  singularities. 

The  elements  of  the  direction  cosine  matrix, 

c-i'yl  ‘-l-  1,2,3.  (2) 

are  nicely  behaved  in  this  sense,  but  they  are  nine  numbers  repre¬ 
senting  three  degrees  of  freedom  and  hence  subject  to  six  con¬ 
straints.  They  do,  however,  possess  one  useful  property  relating 
to  the  structure  of  the  rotations  as  a  group  under  the  operation 
of  "composition"  or  sequential  rotation.  Composition  of  two 
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rotations  is  represented  by  the  matrix  product.  The  rotation 

C3  =  C2  CL  (3) 

is  the  result  of  first  rotating  by  and  then  by  C2.  The  direc¬ 
tion  cosine  matrix  is,  of  course,  similarly  useful  for  mapping 
vectors  in  cartesian  coordinates  directly  (for  further  discussion 
see  Appendix  A) . 

The  Euler  parameters  have  all  the  advantages  of  the  matrix 
product  representation  of  the  group  except  this  point  mapping 
property.  They  are  uniformly  continuous;  that  is,  they  do  not  have 
any  singularities  to  be  contended  with  in  their  calculations  or 
statistics.  Their  primary  disadvantages  are  the  following: 

1.  They  are  obscure, 

2.  They  are  four  parameters  with  one  constraint,  and 

3.  The  representation  is  bivalued. 

The  constraint  is  a  problem  to  be  dealt  with  whenever 
inexact  calculations  of  Euler  parameters  arise:  solutions  to  dif¬ 
ferential  equations,  estimations  based  on  linear  approximations, 
etc.  The  two  representations  of  a  given  rotational  displacement 
differ  only  in  sign,  but  this  is  a  serious  problem  in  numerically 
comparing  two  rotations,  and  it  must  constantly  be  kept  in  mind 
while  developing  numerical  algorithms. 
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2.  DEFINITIONS  AND  FUNDAMENTAL  PROPERTIES 


While  the  motivation  may  not  be  clear,  it  is  convenient  to 
list  here  some  fundamental  terminology.  In  terms  of  the  direction, 
u,  and  magnitude,  6,  of  a  rotational  displacement,  the  Euler  param¬ 
eters  are  given  by 


e0  =  C0S  2 


Bi  =  ui  sin  —  ,  i  =  1,2,3. 


(4) 


For  simplicity,  when  8  appears  without  a  subscript  it  refers  to 
the  column  vector  (8Q,  63,  82,  83).  The  constraint  is 


(5) 


which  is  obvious  from  Eq.  4.  Thus,  the  8^  are  bounded  by  1,  and 

the  bivalued  correspondence  comes  from  the  entrance  of  the  angle 
as  0/2  so  that  a  2ir  rotation  added  to  0  changes  the  sign  of  all 
the 

The  direction  cosine  matrix  is  given  by 


C(8) 


602+B12  ~B22-B32  2(B1B2+B0B3)  2(03^3-0^2) 

2<eiB2-6063>  S0Z-6X2+622-B32  2(B263+60B1> 

_2(B1B3+8082)  2(B2B3-60B1)  B^-B^-b/b-Bj2 


and  it  is  clear  that  C(6)  =  C(-6).  It  is  easy  to  show  that 
(8^,82,83)  is  an  eigenvector  of  C(8)  (with  eigenvalue  1).  This  is 

the  matrix  having  the  successive  transformation  property  described 
by  Eq.  3. 
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A  similar  property  holds  for  the  following  matrix: 


"eo 

-ei 

-g2 

-8 

61 

60 

*3 

-8 

S(8)  = 

e2 

"»3 

eo 

8 

_e3 

62 

“81 

This  matrix  is  introduced  in  Ref.  1  as  having  the  property  (estab¬ 
lished  in  Ref.  2)  that  successive  rotational  displacements  B  and  B' 
result  in  a  rotational  displacement  3"  given  by 


6"  =  S(8')8 


(8) 


in  the  same  sense  that,  as  in  Eq.  3, 


C(8")  =  C(B')  C(8). 


(9) 


Directly,  one  can  show  that  Eq.  8  expands  to 


S(6")  =  S(B')  S (8) ,  (10) 

which  contains  Eq.  8  as  its  first  column.  Thus  it  is  clear  that 
the  Euler  parameters  provide  a  two-valued  representation  of  the 
rotation  group  as  a  group  of  linear  transformations,  one  that 
happens  to  be  a  subgroup  of  the  rotation  group  in  four  dimensions 
(as  represented  by  4  x  4  real  orthogonal  matrices  of  determinant 
1) .  Not  all  orthogonal  4  x  4  matrices  are  of  the  form  of  Eq.  7. 

The  group  identity  is  S(Bq)  =  1^ ,  where  b^  =  (1,0, 0,0). 
And  for  each  8,  the  inverse  of  S(8)  is  its  transpose 

S (8)T  =  S(8) 


Ref.  1.  H.  S.  Morton,  Jr.,  J.  L.  Junkins,  and  J.  N. 
Blanton,  "Analytic  Solutions  for  Euler  Parameters,"  Celestial 
Mech.  Vol.  10  (1974) . 

Ref.  2.  E.  T.  Whittaker,  A  Treatise  on  the  Analytic 
Dynamics  of  Particles,  4th  ed.,  Dover  (1944). 
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where 


K  = 


Then 


R(a)  = 


O’ 

*!•  - 

e2,  -e3) 

the 

matrix 

~-l 

0 

0  0  ~ 

0 

1 

0  0 

0 

0 

1  0 

0 

0 

0  1 

r  any  neR 

4 

> 

KS(a)TK. 

cc„ 

-a. 

-a. 

0 

1 

2 

3 

al 

ao 

"a3 

“2 

a2 

a3 

ao 

~al 

a3 

-a2 

al 

ao 

(11) 


(12) 


(13) 


This  matrix  is  also  given  in  Ref.  1.  It  has  a  property  similar  to 
that  of  S,  Eq.  10:  from 


S(B")  =  S(B')  S(8),  (14) 

we  have 

R(8")  =  KS(6)T  S(6')TK 

=  R(B)R(B' ) ,  (15) 


and  from  column  one,  6"  =  R(8)B*. 


By  direct  evaluation  we  can  establish  the  very  interesting, 
even  amazing,  and  useful  property 

R(a)T  S(8)  =  S(B)  R(cx)T  (16) 

4 

for  any  8  and  any  a  in  R  .  (This  also  implies  that  R(a)  and  S(B) 
commute.)  Next  let  us  observe  that  when  a  =  8,  this  product  takes 
on  the  form 
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1  I  0  0  0 

0  I 

r(6)  =  R(8)  S(8)  =  0  I  c(6)  ,  (17) 

0  i 

where  C(B)  is  the  direction  cosine  matrix  defined  in  Eq.  6. 

This  definition  is  convenient  to  establish  the  law  of  com¬ 
position  stated  in  Eq .  9. 

r(B')E(6)  =  [R(8’)TS(B')]  [R(B)TS(B)1 
=  R(6')T  (R(B)TS(B')1  S(B) 

=  R(B")TS(B") 

=  j:<6")  , 

where  B"  =  S(B')B,  and  Eq.  18  follows  from  Eq.  16,  and  Eq.  19 
follows  from  Eq .  15  (for  R)  and  Eq .  10  (for  S) . 

The  matrix 

T(a,8)  =  R(a) S (8) ,  (20) 

4 

where  a  and  B  are  unit  vectors  in  R  ,  has  six  degrees  of  freedom, 
the  same  as  the  4x4  orthogonal  matrices.  Indeed,  we  can  show  that 
any  orthogonal  4x4  matrix  having  determinant  1  is  of  this  form. 

We  proceed  as  follows:  the  first  column  of  such  a  matrix,  X,  is  a 
unit  vector.  Calling  it  we  form 

S(4)TX  , 

which  then  has  the  form 

1  j  0  0  0" 

o~l - 

o  J  c 
o  l 

T 

where  C  is  orthogonal  and  has  determinant  1.  Thus  S(C)  X  is  of  the 
form  of  Eq.  17: 


(18) 

(19) 
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S(s)TX  =  I(n)  for  some  q.  (21) 

Multiplying  by  S(?)  and  substituting  from  Eq.  17, 

x  =  s(OR(n)Ts(n).  (22) 

We  have  shown  that  S(£)  and  R(n)  commute,  so 

x  =  R(n)Ts[s(c)n] ,  (23) 

and  letting  a=n  (as  defined  in  Eq.  11)  and  B  =  S(c)q,  the  result 
follows.  In  summary,  we  have  shown  that  the  four-dimensional  real 

4 

orthogonal  unimodular  (determinant  +1)  group,  the  rotations  in  R  , 
modulo  (1^,  -1^) ,  is  isomorphic  to  the  direct  product  of  the  groups 

of  S  and  R  matrices,  modulo  (1^,  -1^). 
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3.  QUATERNIONS  AND  OTHER  REPRESENTATIONS  OF 
THE  ROTATION  GROUP 

We  have  dealt  up  to  this  point  with  several  representations 
of  the  rotation  group:  the  direction-cosine  matrices,  the  S  and  R 
matrix  representations  that  we  associated  with  the  term  "Euler 
parameters,"  and  others  were  mentioned.  Some  authors  prefer  to  call 
the  Euler  parameters  quaternions,  and  indeed  the  representations  are 
isomorphic.  Historically,  the  parameters  were  first  developed  by 
Euler  and  Rodrigues  and  were  published  in  1776.  They  have  variously 
been  called  "the  homogeneous  parameters  of  Euler,"  "Euler's  sym¬ 
metric  parameters,"  and  "Euler-Rodrigues  parameters."  It  was  in 
1843  that  Hamilton  introduced  quaternions  and  the  properties  of  this 
division  algebra,  not  subject  to  the  length  constraint. 

Figure  1  illustrates  the  diversity  of  terminology  in  the 
standard  texts  and  technical  literature  that  has  evolved  over  the 
last  200  years.  Keep  in  mind  that  all  of  this  addresses  one  pri¬ 
mary  mathematical  object:  the  group  of  rotations  in  3-space. 
Representations  linked  by  double  arrows  are  identical  or  nearly 
identical.  Numbers  identify  references  at  the  end  of  this  report. 
Arrows  without  numbers  can  mostly  be  found  in  Chapter  4  of  Ref.  5. 
Other  direct  transformations  beyond  those  shown  by  the  arrows  cer¬ 
tainly  exist  or  can  be  derived,  for  example,  the  calculation  of 
z-x-z  Euler  angles  from  direction  cosines.  To  make  matters  worse 
or  at  least  more  confusing,  there  is  no  universal  standard  nota¬ 
tion  within  individual  representations.  Euler  angles,  for  example, 
may  be  defined  in  many  different  ways. 

By  the  term  "principal  axis  form"  is  meant  a  specification 
of  the  axis  of  a  rotation  and  an  angle.  This  can  be  a  unit  vector 
and  a  scalar  for  the  angle,  or  the  product  of  the  angle  and  the 
unit  vector.  The  Rodrigues  parameters  are  not  discussed  in  any  of 
the  references,  but  they  are  essentially  like  the  principal  axis 
form:  the  unit  vector  is  multiplied  by  the  tangent  of  one-half  of 
the  rotation. 

A  representation  of  a  group  must  model  the  group  operation. 
The  ease  with  which  this  is  done  provides  one  distinction  between 
representations.  The  direction  cosines  and  Euler  parameters  take 
the  form  of  real  matrices  and  combine  according  to  the  ordinary 

Ref.  3.  G.  Birkhoff  and  S.  MacLane,  A  Survey  of  Modern 
Algebra,  MacMillan  (1953) . 

Ref.  4.  H.  C.  Corben  and  P.  Stehle,  Classical  Mechanics, 
2nd  ed.,  Wiley  (1950). 

Ref.  5.  H.  Goldstein,  Classical  Mechanics,  Addison- 
Wesley  (1950) . 
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Fig.  1  Representations  of  the  group  of  rotations. 
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matrix  multiplication.  The  Pauli  spin  and  Caley-Klein  representa¬ 
tions  use  complex  matrices  of  lower  order  and  are  basically  equiva¬ 
lent  in  computational  complexity  to  quaternions  and  Euler  param¬ 
eters.  Principal  axis  and  Rodrigues  notation  use  vector  cross  and 
dot  products,  and  Euler  angles  require  complicated  trigonometric 
calculations  plus  logical  branching  to  avoid  problems  at  singular¬ 
ities. 
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4.  KINEMATICS  AND  DYNAMICS  OF 
ROTATIONAL  MOTION 


From  Eq.  10  come  the  interesting  properties  of  the  Euler 
parameters.  Differentiating  Eq.  4,  we  obtain 


(24) 


This  form  can  be  simplified  a  bit  by  noting  that  when  6  =  0  we  have 


(25) 


where  =9^,  which  are  the  three  components  of  the  instantaneous 
spin  vector.  Next  we  can  differentiate  Eq.  10  to  obtain: 

S(8”)  =  S(8’)S(8)  +  S(6')S(B).  (26) 


Now  let  8  be  a  constant,  selected  to  be  equal  to  8"  at  some  time,  t^. 
Then  8  =  0  and  at  tg,  8"  =  8  so  that  at  tg  we  have  8'  »>  bg,  the  iden¬ 
tity.  Thus  the  above  form  applies  to  6*  and  Eq.  26  becomes 


S[8"(t0)]  =  S(«/2)Sl8"(t0)J  =  y  S(u>)St8"(t0)].  (27) 
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The  first  colvui  of  Eq.  27  becomes 
B"(t0)  -  ^S[u(t0)38"(t0), 
which  holds  for  any  choice  of  tgj  thus  at  t: 

0(t)  =  yS[w(t)]6(t),  (28) 

where  w(t)  =  0 

ux(t) 
u2(t) 

_“3(t) 

and  the  last  three  components  are  the  instantaneous  spin  vector  in 
body-fixed  coordinates.  Thus 

S[0(t)]  =  |s[u(t)]S[B(t)]  .  (29) 

The  reason  u  appears  in  body-fixed  coordinates  is  that  the 
spin  motion  follows  the  displacement  8.  This  is  due  to  the  choice 
in  Eq.  27  to  hold  6  constant.  Had  we  held  0'  constant  and  let  0 
vary,  the  spin  in  inertial  coordinates  would  appear. 

Also  we  see  that  the  form  of  Eq.  10  is  that  of  the  state 


transition  matrix: 

8(0  =  *(t,t0)B(t0)  (30) 

in  Ref.  1.  Let  us  define,  for  any  8(t),  a  matrix 

♦(t,t0)  =  s[0(t)]s[e(to)]T  ,  (3i) 

which  is  a  rotational  motion,  the  inverse  of  0(tg)  followed  by  8(t). 
Call  its  first  column  a(t),  so  that 

*(t,t0)  =  Sfa(t)]  (32) 

is  defined  properly,  and  observe  from 

S[0(t)]  -  S[a(t) ]S[0(tQ) ]  (33) 


20  - 


THE  JOHNS  HOPKINS  UNIVERSITY 

APPLIED  PHYSICS  LABORATORY 

LAUREL.  Maryland 


that  Eq.  30  does  indeed  hold  for  this  choice  (Eq.  31)  of  (Kt.tg). 
Differentiating  Eq.  33  and  applying  Eq.  28  we  obtain 

yS[u(t)]S[B(t)]  =  *(t,t0)S[B(t0)]  , 
or 

*<t,to)S[0(to)]  =  S[u(t)]*(t,t0)S[B(t0>]  , 

and  thus 

$(t,tQ)  =  yS[w(t) ]$(t,t0)  ,  (34) 

as  the  differential  equation  for  the  transition  matrix.  All  of  the 
preceding  is  kinematic  and  is  true  for  any  history  of  w(t). 

We  can  use  Eq.  17  in  developing  the  dynamics  of  rigid  body 
motion.  We  need  the  derivative  of  C(B)  with  respect  to  time  in 
terms  of  the  body-fixed  spin  components,  u. 

|ti(B)  =  i(B)  =  R(B)XS(B)  +  R(B)TS(B)  .  (35) 

Because  R  and  S  are  linear  in  B  (E  is  quadratic), 

E(8)  =  R(B)TS(B)  +  R(B)TS(B)  .  (36) 

The  form  of  interest  is 


Z(B)E(6)T  =  R(B)TS(B)S(B)TR(B)  +  R(B)TS(B)S(B)TR(B) 


=  R(B)TR(B)  +  S(B)R(B)TR(6)S(B)T 

(37) 

=  R(B)TR(B)  +  S(B)  S(B)T  . 

(38) 

We  know 

tion  of 

•  T 

S(B)  from  Eq.  29  and  using  Eq.  12  we  find  R(6)  . 
these  yields  yS(«)  for  each  term  in  Eq.  38,  so 

Evalua- 

i(B)r(B)T  =  S(w)  , 

(39) 

where  w 

=  (0,  0)^,  oo2>  w3)  in  this  notation. 
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Is 


The  fundamental  physical  law  of  rotation  of  rigid  bodies 


d 

dt 


(40) 


where  t  is  the  torque  applied  to  the  body  and  L  is  the  angular 
momentum.  The  subscript  "i"  reminds  us  that  this  is  true  in  iner¬ 
tial  reference  frames.  We  can  define  the  inertia  tensor 


J  = 


/ 

body 


(rTrI-rrT)  pdV 


(41) 


in  which  dV  is  differential  volume,  p  is  density,  I  is  the  3x3 
identity,  and 


r  = 


is  a  point  in  the  body. 


T  T 

The  matrix  (r  rl-rr  )  is  then 


(42) 


(43) 


which  is  valid  only  in  the  specific  coordinate  frame  in  which  r  is 
defined.  Thus  a  rotating  body  will  have  an  inertia  tensor  that 
may  be  time-varying  in  any  coordinates  other  than  body-fixed. 

Angular  momentum  is  conveniently  separated  as  the  product 

L  *=  Jio,  (44) 


and 


L  =  jo)  +  Ju.  (45) 

As  an  alternative  to  writing  J,  we  can  transform  Eq.  40  into  body- 
fixed  coordinates  using 

Jt  -  C(8)T  Jfc  C(8)  ;  (46) 
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«0±  =  C(8)T  ub  ;  (47) 

t±  =  C(8)T  xb  .  (48) 

•  • 

We  must  substitute  for  J  and  w  in  Eq.  45 

jj.  =  C(8)TJb  C(8)  +  C(6)T  Jb  C(B)  ;  (49) 

U).  =  C(8)T  ub  +  C(8)Tii>b  .  (50) 

This  produces,  in  place  of  Eq.  45 


Li  =  C(8)TJbwb  +  C(6)T  Jfa  C (6)  C(8)T 

+  C(8)T  Jb  C(8)  C(6)T  ub  +  C(6)T  Jb  wb  .  (51) 

Equating  to  torque  (Eq.  40)  and  substituting  with  Eq.  48 
to  involve  body-fixed  representation  exclusively,  we  obtain 

L±  =  Tt  =  C(8)T  Tb  (52) 

and 

C(B)  La  =  rfe  .  (53) 

Multiplying  Eq.  51  on  the  left  by  C(8)  simplifies  it  somewhat  and 

•  •  T 

and  makes  all  occurrences  of  C  be  of  the  form  CC  or  the  transpose 
of  this.  From  Eq.  39  we  can  extract  the  lower  right  3  x  3  to  leave 
C  in  place  of  E,  and  write  (calling  this  submatrix  of  S(w)  by  the 
same  name) , 

Tb  =  s(“b)T  Vb  +  Jb  s(aib)wb  +  Jb  s(wb)T  “b 

+  Jb  “b  •  (54) 

Note  that  wb  is  used  in  S  because  we  have  chosen  in  Eq.  46  through 
Eq.  48  the  transformation  C(B)  to  be  from  inertial  to  body-fixed. 
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Next,  observe  that  Eq.  54  simplifies  by  skew- symmetry  of 

S(w,  )  to  the  familiar 

□ 

Vb  =  S(V  Jb“b  +  Tb  *  (55) 

This  is  the  body-fixed  form  of  Eq.  40  with  j  terms  in  more  tract¬ 
able  form. 

Equations  55  and  28  or  29  form  a  set  of  differential  equa¬ 
tions  for  the  motion  of  a  rigid  body.  In  this  context  we  need  to 
indicate  the  dependence  of  t  on  the  other  variables,  specifically 
8,  a),  and  time: 

Tb  =  Tb  *  (56) 

Nowhere  else  in  Eq.  55  does  8  appear,  and  thus  a  distinction  can 
be  made  between  two  classes  of  problems:  those  in  which  t  depends 
on  8  and  those  in  which  it  does  not.  For  numerical  purposes,  it  is 

a:  that  is  of  interest,  and  it  always  depends  on  w  except  in  trivial 
problems . 

Torques  that  do  not  depend  on  6  are  those  due  to  on-board 
rockets,  momentum  wheels,  fluids,  nonrigidity,  or  other  mass-move¬ 
ment  effects.  Those  torques  derived  from  interaction  with  external 
masses  or  fields  (drag,  radiation  pressure,  magnetic,  gravity  gra¬ 
dient)  do  depend  on  8.  The  simplification  in  the  former  case  al¬ 
lows  solution  numerically  to  the  differential  equation  (Eq.  34) 
for  the  entire  state  transition  matrix  without  knowledge  of  the  ini¬ 
tial  attitude,  B(tp).  This  problem  is  a  first-order  differential 

equation  (nonlinear)  in  three  variables;  Eq.  34  for  8(t)  is  reduced 
a  quadrature. 

There  are  important  implications  here  for  satellite  atti¬ 
tude  determination  in  the  sense  that  the  transition  matrix,  4>,  for 
the  8  variables  is  itself  an  S-matrix  (Eq.  32)  and  can  be  repre¬ 
sented  as  S[a(t)]  for  some  time-varying  Euler  parameters,  a(t). 

Then  the  matrix  C[8(tp)]  can  be  propagated  as  well: 

C[8(t) ]  =  C[a(t) ]  C[8(t0) ]  ,  (57) 

and  this  inverted  to 

C[8(t0)]  =  C[a(t) ]T  C [ 8 ( t ) ] 


(58) 
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to  translate  back  to  t^  any  future  information  about  the  satellite 

attitude  in  the  form  of  observation  of  directions  in  space.  Most 
attitude  sensors  (magnetometers,  solar  detectors,  star  cameras) 
that  are  not  rate  sensors  are  of  this  type. 

The  implication  of  Eq .  58  specifically  is  this:  in  the 
absence  of  8-dependent  torques,  a  knowledge  of  the  initial  spin, 
u(tp),  reduces  the  attitude  determination  problem  to: 

1.  A  first-order  differential  equation  in  to  only, 

2.  A  quadrature  to  yield  C[ot(t)],  and 

3.  A  least-squares  estimation  of  8(tg). 

That  is,  the  time  dependence  is  removed  from  the  estimation  process; 
it  is  as  if  all  measurements  are  taken  simultaneously  at  t^.  One 

can  even  propagate  the  covariance  effects  of  errors  in  the  knowledge 
of  to  back  from  the  time  of  each  measurement  to  derive  appropriate 
weights  for  the  estimation  of  S  at  tQ. 

Numerical  solutions  of  Eqs.  28  and  56 

8  =  \  S  (to)  8  , 

to  =  J  ^  S(to)  Jto  +  J  ^  x(6,u),t)  >  (59) 

may  be  improved  considerably  in  computation  time  by  the  following 
method  in  cases  when  to  varies  much  more  slowly  than  8.  In  general 
this  will  include  only  problems  in  the  class  having  torques  inde¬ 
pendent  of  8.  In  this  case  there  is  much  to  be  gained  by  replacing 


the  approximation 

A8  =  8  At 

(60) 

by  yS  (to)  At 

B(t+At)  =  e  6(t)  , 

(61) 

where  8  and  to  denote  values  at  time  t,  which  we  fix  at  the  begin¬ 
ning  of  the  interval  of  length  At.  The  matrix  exponential  serves 
the  role  of  the  transition  matrix  when  to  is  constant: 

\  S(to) At 

$ (t+At , t)  =  e  (62) 

as  in  Eq.  30  and  satisfies  the  differential  equation  (Eq.  34)  when 
the  derivative  is  taken  with  respect  to  At  holding  t  (and  to) 
constant. 
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It  is  not  necessary  to  use  the  series  expansion  for  the  matrix  ex¬ 
ponential  or  even  to  retain  the  exponential  notation,  for  we  may 
refer  back  to  the  fundamental  definitions  (Eqs.  4  and  7)  and  write 


4>(t4At,t)  =  S(cos 


u  sin 


to  |  At  v 
2  '  ’ 


(63) 


where  u  is  the  unit  vector  as  given  in  Eq.  4  with  an  appropriately 
chosen  sign.  Appendix  B  shows  a  different  approach  to  the  develop¬ 
ment  of  Eqs.  62  and  63. 


A  Runge-Kutta  fourth-order  algorithm  based  on  Eq.  63  was 
developed  by  straightforward  replacement  of  all  B(t  +  At)  =  8(t)  + 


6(t)At  or  similar  steps  by  the  rotational  equivalent.  Some  atten¬ 
tion  is  also  required,  however,  to  the  order  of  matrix  products, 
since  these  are  no  longer  commutative  operations.  Let  h  be  the 
time  step  and  define 


0(u>,At)  =  S 


cos 


1  fa)  1  At 
2 


(64) 


Thus  9(fa), At)  represents  the  rotation  generated  over  an  interval  At 
by  the  constant  spin  vector,  u>.  The  Runge-Kutta  algorithm  is  as 
follows: 


\  =  0[u(ti),  h/2] 

(65) 

kl  =  f  e(ci) »  “(ti) ,  t± ] 

(66) 

6’  =  e ( t ± ) 

(67) 

w'  =  k^  +  w^) 

(68) 

X2  +  0(fa)' ,  h/2) 

(69) 

k2  =  J  T(6' »  +  h/2) 

(70) 

6*  =  X2  6(t±) 

(71) 

u*  «  k2  +  ^(t^) 

(72) 
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A3  =  0(w*.  h)  (73) 

k3  =  h  t(6*,  w*.  t±  +  h/2)  (74) 

@+  =  A3  6(t.)  (75) 

w"t"  =  k3+w(ti)  (76) 

k4  =  |  T(e+,  u>+,  t.  +  h)  (77) 

3(ti+1)  =  0(c+,  h/6)  0(^y^-,  h/3)  e[U(t.),h/6]  f$(t.) 

(78) 

u(ti+l)  =  “(ti}  +  (ki  +  2k2  +  k3  +  k4)/3  .  (79) 


When  the  torque,  x,  does  not  depend  on  6,  it  is  possible  to 
to  modify  this  algorithm  to  integrate  0(t)tp)  from  the  differential 

equation  (Eq.  34)  and  the  initial  condition  $(t0,t0)  =  1^,  the  4x4 

identity  matrix.  The  intermediate  calculations  of  6  would  be  un¬ 
necessary  and  |t(ti+3,tQ)  would  be  produced  in  Eq .  78  by  replacing 

6(ti)  by  $(ti,  tQ). 

As  a  numerical  example,  the  torque-free  motion  of  a  body 
whose  inertia  moments  were  =  43.1,  J33  =  40.6,  and  J33  =  44.3 

was  integrated  from  an  initial  6(tp)  =  (1,  0,  0,  0),  and  w(tp)  =  10 

rpm,  offset  by  0.1  radian  (about  6°)  from  the  z-axis  in  the  +x 
direction. 

Figure  2  shows  the  spin  history  for  500  sec,  about  83  revo¬ 
lutions.  Precession  has  carried  the  spin  axis  through  four  trips 
around  the  z-axis.  Using  the  Runge-Kutta  scheme  presented  here, 
one  can  select  a  step  size  h  based  on  this  lower  rate  of  u>  vari¬ 
ation  rather  than  the  high  rate  of  6  periodicity.  Variations  of  6  j 

are  shown  over  50  sec  in  Fig.  3,  Table  1  shows  roughly  the  order  i 

of  magnitude  of  accumulated  errors  in  6  and  in  u  at  the  end  of  500  j 

sec.  These  values  were  obtained  by  comparison  of  results  with  values  J 

obtained  by  integrating  at  2  sec  steps.  j 

I 
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Table  I 

Krrors  in  w  and  d  from  numerical  sources. 


Step 

size 

(sec) 

5 

10 

20 

Au 

0.0001 

0.001 

0.015 

A6 

0.008 

0.02 

0.10 

The  errors  in  6  fall  within  the  range  from  ^  to  — , 

where  they  are  expected  to  be  if  they  are  due  to  errors  in  u.  Thus, 
one  would  conclude  in  this  case  that  the  8  integration  is  "more 
linear"  than  the  w  integration.  The  opposite  would  be  true  in  the 
case  of  large  but  constant  (in  body  coordinates)  torques. 
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5.  A  PROBABILITY  FUNCTION  ON  THK  CROUP 


The  size  or  "amount"  of  rotation  indicated  by  a  vector,  6, 
of  Euler  parameters  can  be  measured  directly  by  calculating  the 
angle  8  in  the  definition  given  by  Eq.  4.  A  comparison  of  the 
"difference"  between  two  orientations,  8'  and  B,  can  be  made  by 
using  the  group  operation  and  inverse  as  a  means  of  "subtraction." 
On  writing 

8"  =  S(S’)T  8  ,  (80) 


we  could  calculate  the  angle  9"  that  corresponds  to  8";  however, 
there  are  equivalent  measures  of  the  amount  of  rotation  that  are 
simpler  to  calculate.  For  example,  we  can  define  a  metric  on  the 
rotation  group  as  follows: 

P (6, 8' )  =  | sin  8"/2|  ,  (81) 


where  0"  is  as  defined  above.  The  reason  this  is  simpler  is  that 
it  eliminates  trigonometric  calculations  by  simplifying  to 


P (6, 8 ' )  = 


£ 


-  (bV)2 


(82) 


This  can  be  seen  by  observing  that  the  first  element  of  Eq.  80, 

T  Q  ** 

8q",  is  8  6',  which  is  also  cos  — ^  .  Thus 


(cos  -^)2  =  (SV)2 


(83) 


and  Eq.  82  follows  directly.  Note  that 

P (8,-8)  =  0  ,  (84) 

which  means  that  p  is  a  metric  on  the  rotation  group  and  not  on  the 
4 

unit  sphere  in  R  because  it  does  indeed  identify  antipodal  points. 
We  have  not  shown  that  p  satisfies  the  definition  of  a  metric,  in 
particular  the  triangle  inequality.  This  calculation  is  tedious 
but  not  difficult. 


An  important  observation  is  the  invariance  of  p  under  rota¬ 
tion, 


p [S(a)8,  S (a) 8 '  ]  =  p(6,8’) 


(85) 
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T 

which  follows  directly  from  the  invariance  of  6  B'  under  such  an 
operation.  It  is  also  interesting  to  note  that  the  function  p  can 
be  written  as 


4 

and  extended  to  all  of  R  ,  where  the  projection  back  to  the  unit 
sphere  is  incorporated  in  p  directly.  Thus,  one  should  beware  of 
numerical  errors  in  the  length  of  f?  when  applying  p  in  the  form 
Eq.  82,  which  is  not  corrected  for  length. 

This  metric  notion  of  the  amount  of  rotation  leads  directly 
to  the  concept  of  the  scalar  magnitude  of  a  rotation  defined  as 

||6l|  =  p(B,b0)  ,  (86) 

so  that  (|bg!|  =0,  || e|f  sl  for  any  S  and,  from  Eqs.  81  and  4, 

||s!|  -  |  sin  0/2)  .  (87) 

We  can  speak  of  a  differential  rotation 


ee 


( 

I 


.  de 

sin  — j 


i  =  1,2,3 


(88) 


as  in  Eq.  4,  which  has  a  differential  scalar  magnitude  although 

4 

it  is  a  unit  vector  in  R  .  The  first  component  is  large  and  is 
not  deQ  in  the  usual  linear  sense;  for  this  reason  we  choose  the 

Greek  6  rather  than  d,  emphasizing  this  distinction. 
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The  calculus  of  these  differential  rotations  and  of  distri¬ 
butions  or  functions  defined  on  the  rotation  group  can  be  extracted 
from  the  more  general  and  abstract  methods  presented  in  Refs.  6 
through  10;  however,  in  an  effort  to  gain  clarity  and  insight,  we 
choose  to  develop  methods  in  the  specific  context  not  only  of  the 
rotation  group  but  of  the  particular  representations  of  it  that  are 
most  familiar,  such  as  Euler  angles,  Euler  parameters,  etc. 


Let  us  denote  by  G  the  rotations  as  an  abstract  group,  inde¬ 
pendent  of  any  representation  by  parameters.  Suppose  f  is  a  real¬ 
valued  function  defined  on  G.  We  will  find  it  useful  to  be  able  to 
integrate  f  over  all  or  part  of  G,  such  as  one  integrates  probabil¬ 
ity  densities  to  obtain  probabilities.  But  the  symbols 


f (g)dg 


(89) 


G 


are  not  well  defined  without  some  definition  of  what  is  meant  by 
dg. 

We  need  a  "measure"  on  G,  and  the  symmetry  or  uniformity  of 
G  ought  to  imply  that  for  any  fixed  rotation  gQeG, 

/f(gg0)dg  =  (g)dg  ,  (90) 

G  G 


that  is,  that  this  measure  is  uniform  over  G.  Take,  now,  as  an 
example  a  representation  of  G,  the  z-x-z  Euler  angles,  which  we 
denote  by  (<Ji,0,<|>).  Write  g  as  a  function  of  these  so  that  f(g)  = 
f(<j>,0,i{/)  and 


f  f(g)dg 


2tt,tt,2it 


f(M,*)K*,e.*)d(M0di|» 

0,0,0 


(91) 


Ref.  6.  I.  M.  Gel' f and  and  Z.  Ya.  Sapiro,  "Representations 
of  the  Group  of  Rotations  in  Three-Dimensional  Space  and  Their  Ap¬ 
plications,"  AMS  Trans . ,  Series  2,  Vol.  2  (1956). 

Ref.  7.  F.  D.  Murnaghan,  The  Theory  of  Group  Representa¬ 
tions,  Dover  (1963)  (alco  pub.  by  Johns  Hopkins  Press,  1938). 

Ref.  8.  U.  Grenander,  Probabilities  on  Algebraic  Struc¬ 
tures,  Wiley  (1963) . 

Ref.  9.  C.  Chevalley,  Theory  of  Lie  Groups,  Vol.  I, 
Princeton  (1946) . 

Ref.  10.  H.  Flanders,  Differential  Forms,  Academic  Press 

(1968). 
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If  we  can  find  a  "weight  function," 
and  normalize  it  so  that 


8tt 


I, 


with  the  property,  Eq.  90, 


(92) 


where 


dg  ■»  I (4> ,  e , <P) d<t>d6d(jj  ,  (93) 

then  Eqs.  91  and  93  can  become  a  definition  of  dg  that  makes  Eq.  89 

2 

well  defined.  The  choice  of  8tt  will  be  clear  later. 

Intuitively,  we  see  that  for  the  differential  rotations 

(d<(>,d0,diji)  to  be  orthogonal,  we  must  have  0  =  y,  for  which  d<j>  is 

about  z,  d0  about  x,  and  dif/  about  y.  The  alignment  of  the  if>  and  ip 
angles  when  0  is  small  means  that  the  l($,0,if»)  value  ought  to  be 
small  too.  As  it  turns  out,  I(i|i,0,iJ>)  =  sin  0  works.  In  this  ex¬ 
ample,  the  key  feature  of  tnis  weighting  function  is  that  it  repre¬ 
sents  the  nonorthogonality  of  the  three  differential  rotations. 

The  invariant  (under  rotation)  quantities  of  interest  are  the  "dif¬ 
ferential  volume,"  represented  by  differential  rotations  in  three 
degrees  of  freedom,  and  the  "relative  orthogonality"  of  these  axes 
or  "angles"  between  them.  The  function  I((£,8,<Ji)  =  sin  0  represents 
the  volume  of  a  parallelepiped  whose  edges  are  unit  vectors  in  the 
directions  of  these  three  Euler  rotations  as  viewed  in  a  fixed  (not 
the  running  z-x-z)  coordinate  frame.  It  corrects  the  quantity 
di(id0di>  to  the  differential  volume,  dg,  that  is  needed  as  an  invari¬ 
ant  under  rotation  so  that  the  integral  satisfies  Eq.  90. 

We  can  do  the  same  thing  using  the  Euler  parameters.  Since 
there  are  four  of  these  and  only  three  degrees  of  freedom  in  the 
rotation  group,  the  differential  volume  will  be  a  "three-form"  in 

4 

the  exterior  algebra  over  R  ,  an  elaborate  way  of  saying  that  the 
previous 


dg  *■  I($, 0,i|Od<j>d0dij> 
will  be  replaced  by 

dg  -  I0(8)d81d62d83  -  I1(6)d0od62de3 

+  I2(8)d80d61d63  -  I^dBgdB^ 
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and  1(8)  (a  vector)  is  to  be  chosen  to  make  dg  an  invariant  under 
rotations.  For  notational  simplicity  Eq.  94  will  be  written 

dg  =  I(8)T  *dB  ,  (95) 

where  *d3  denotes  a  vector  (*  is  the  Hodge  star  operator,  Ref.  10), 


*d8  = 


'  de1de2d63 

"dW63 

d80d81d83 

-d60d81d82 


(96) 


Next,  from  Eq.  88  we  see  that  at  8  =  b^  =  (1,0,0, 0)  we  ought 
to  have  dg  =  8  dB^dB^B.^.  Thus  we  conclude  we  want  1(1, 0,0,0)  = 

(8, 0,0,0).  The  requirement  of  invariance  is,  in  our  current  notation, 

KB')1  *d8'  =  I ( 6) T  *dB  (97) 

for  any  choice  of  8  and  8'.  Fortunately,  we  know  how  Euler  param¬ 
eters  are  related;  there  exists  an  a  such  that 

8'  =  S(ot)8  .  (98) 


From  this  it  follows  that,  considering  a  a  constant, 

dB'  =  S(a)d8  (99) 

and 


*dB’  =  S(c»)*d3  ,  (100) 

which  can  be  argued  from  the  framework  of  exterior  algebra  in  which 
and  "d"  are  linear  operators,  or  by  accepting  Eq.  99  on  the 
basis  of  "d"  being  linear  and  calculating  Eq.  100  from  the  defini¬ 
tion  of  *d3. 

In  any  event,  multiplying  Eq.  100  on  the  left  by  the  1x4 
matrix  I(B')T, 

I(B')T  *dB'  =  I(B')T  S(a)  *dB  ,  (101) 

and  thus 

I (8)T  =  I(8’)T  S(a)  ,  (102) 
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I(B')  =  S (a)  1(8)  ,  (103) 

which  holds  for  any  8,  specifically  8  =  (1,0, 0,0).  So  we  obtain  by 
substituting  this  8,  and  8'  =  a, 

1(a)  =  8a  .  (104) 

Thus  the  invariant  differential  volume  element  is 


dg  =  8(80d81d82d83  - 

+  82d80d31d83  -  B^BgdB^B^ 


(105) 


This  differential  form  is  defined  on  the  surface  of  the  unit  ball 

4 

in  R  ,  so  we  can  apply  Stokes'  Theorem  and  actually  evaluate  the 
integral  to  obtain 


/dg  -  32  / 

]e|=i 


d80d81d6;;!d83  =  16ir 


8  <  1 


(106) 


This  integral  is  twice  as  large  as  Eq .  92  because  of  the  bivalued 
representation  of  the  rotation  group;  that  is 


(107) 


as  in  Eq.  92.  Note  that  the  topological  and  geometric  properties 
result  in  expressions  symmetric  in  the  four  components  of  8  (Eqs. 
105  and  82);  the  uniqueness  of  8q  as  in  Eq.  4  appears  in  expres¬ 
sions  relating  to  the  group  structure  and  in  particular  the  group 
identity. 


Having  thus  demonstrated  an  invariant  volume  element,  dg, 
we  can  think  about  probability  densities  on  G,  those  functions,  f, 
normalized  to  have  their  integral  Eq.  89  over  the  group  be  unity. 
For  example,  the  constant 

f(g)  -  -K 

8* 


is  one  such  and  represents  the  "equal  likelihood"  density  function, 
or  the  "uniform  density  function."  The  finite  volume  of  the  group 
and  the  existence  of  a  uniform  density  function  do  not  parallel  the 
structure  of  the  real  line  and  the  theory  of  probability  distribu- 
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tions  of  real  variables.  Another  major  distinction  is  that  the 
properties  of  linear  spaces  are  enjoyed  by  the  real  line  and  not  by 
the  rotation  group.  Both,  however,  have  a  group  structure  that  al¬ 
lows  the  fundamental  operation  of  convolution  to  be  defined.  A 
theory  of  Fourier  transforms  can  be  developed  for  G  as  it  can  for 
the  real  line  (Refs.  7  and  8). 


We  are  going  to  investigate  an  analog  of  the  multivariate 
normal  density  on  Euclidean  space,  the  probability  density  function 


PD(*) 


k(D)  e 


1  T 
j  6  DB 


(108) 


where  D  is  a  real  4x4  symmetric  matrix,  and  k(D)  is  a  normaliza¬ 
tion  constant  defined  by 


1  = 


/ 


PD(6)dg  =  j  k(D) 


/ 


(109) 


There  are  two  primary  motivations  for  dealing  with  this  choice: 
first,  given  two  such  functions,  the  product  is  another  such  after 
normalization  and  admits  a  maximum-likelihood  estimation  analogous 
to  that  for  the  multivariate  normal  probability  functions;  and 
second,  densities  of  this  form  represent  well  the  information  ob¬ 
tained  from  single  measurements  of  commonly  used  attitude  determin¬ 
ation  instruments  on  spacecraft. 


A  number  of  facts  are  immediately  obvious.  The  choice  of  D 
to  represent  the  density  function  is  not  unique.  Let  D'  =  D  +  XI, 

X  a  constant,  then 

BTD'B  =  X  +  eTDB  ,  (110) 


and 


PD, (6)  =  k(D') 


— X / 2  “  2  ^ 

e  e 


(111) 


Therefore,  pD (3)  and  pD,(B)  differ  at  most  by  a  constant,  but  by 
Eq.  109  they  must  be  equal: 

pD<8)  =  PD'<6) 

and  (112) 

k(D')  =  k(D)  e  ' 
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We  also  can  observe  that  the  uniform  distribution  is  a 
member  of  this  class  of  functions.  Letting  D  =  XI,  the  dependence 
on  8  is  eliminated,  for  example,  with  X  =  1: 

-  eTe  -  - 

Pl(6)  =  k(I)  e  "  2  =  k(I)  e  2  ,  (113) 

and  thus 

e1/2 

k(I)  =  S—  .  (114) 

8ir 

What  we  have  is  a  one-to-one  correspondence  between  these  density 
functions  and  the  equivalence  classes  of  real  symmetric  4x4 
matrices  defined  by  the  equivalence  relation  that  two  matrices  are 
equivalent  if  their  difference  is  a  constant  times  the  identity 
matrix.  This  is  the  reason  we  have  not  mentioned,  for  example, 
positive  definiteness;  it  is  not  a  class  property,  although  each 
class  clearly  has  a  positive  definite  representative. 

We  cannot  tie  estimation  to  the  "expectation"  of  8  (nor 
even  define  it)  because  G  is  not  a  linear  space,  so  we  must  appeal 
to  maximum  likelihood  for  an  estimation  of  8.  If  we  were  dealing 
with  a  probability  over  the  full  four-dimensional  space,  we  could 
solve  for  the  peak  value  of  p^  (the  most  likely  8)  by  differenti¬ 
ating  with  respect  to  8  and  setting  this  equal  to  (0,0, 0,0). 


d  p  (8)  -  \  6TD8 

— ft -  -  -k(D)  e  L  (8  D)  .  (115) 

T 

But  this  leads  only  to  8  D  =  0  and  8=0,  which  is  not  very  useful. 
Incorporating  the  constraint  |sj  =1,  however,  we  can  find  the  de¬ 
sired  so]  ition  by  determining  points  for  which  the  gradient  of  p^ 

is  parallel  to  8.  From  among  these  we  select  the  one  for  which  pD 

is  largest.  When 

d  PD(8) 

~ dB  =  C  6  ,  (116) 


with  c  constant,  the  directional  derivatives  perpendicular  to  8  are 
0,  these  being  the  only  directions  that  preserve  |8|  =1. 
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This  equation  is  written,  from  Eqs.  115  and  116,  as 

-  Pd(B)  (6TD)  =  C8T  .  (117) 

Solutions  of  this  equation  are  the  eigenvectors  of  D,  for  which 
there  are  four  real  eigenvalues,  X^.X^.X^.X^,  Let  a  set  of 

orthogonal  eigenvectors  corresponding  to  these  be  a^,.,.,ot^  and 

suppose  the  X^  are  denumerated  in  ascending  order.  We  now  have 

-  K 

PD(ai)  =  k(D)  e  ,  (H8) 


so  Pp  is  largest  for  the  smallest  of  the  namely  X^,  although  ot^ 
may  not  be  unique,  as  is  the  case  for  p^,  for  example. 

Another  important  property  is  that  a  probability  density 
function  may  be  "rotated."  If  we  form 

-  j  BTS(cOTDS(a)8 

PD[S(a)8]  =  k(D)  e  44  ,  (119) 

we  obtain  Pg,(a),  where 

D'  =  S(a)T  DS (a)  (120) 

and  by  the  invariance  of  the  normalizing  integral, 

k(D’)  =  k(D)  .  (121) 

That  D’  is  symmetric  is  obvious;  it  clearly  has  the  same  eigen- 

T 

values  as  D  and  eigenvectors  S (ct )  ou  corresponding  to  them.  For 
example,  setting  a  =  makes  bg  [bg  =(1,0, 0,0)]  the  most  likely 
orientation. 

As  was  previously  mentioned,  the  product  of  pg  and  Pg,  is  a 
function,  Pg„,  of  the  same  type: 

-  •Wr(D+D,)f3  .  (122) 

Pg,(8)  Pg(p)  -  k(D')  k(D)  e 


On  writing 


D"  -  D+D’ 


-  39  - 


(123) 


THE  JOHNS  HOPKINS  UNIVERSITY 

APPLIED  PHYSICS  LABORATORY 

LAUREL  MARYLARO 


we  find  pD„  is  normalized  according  to 


PD-  (6) 


k(D+D') 
k(D)  k(D') 


V(6) 


pD(6) 


(124) 


There  is  a  similarity  here  to  the  process  of  combining  two  linear 
least-squares  estimates  based  on  independent  sets  of  observations. 
These  "D"  matrices  play  the  role  of  the  "information  matrix"  or  the 
inverse  of  the  "error  covariance  estimate."  It  is  the  differences 
rather  than  the  similarities  that  are  of  interest.  As  we  have 
shown,  the  "zero  information"  matrix  is  not  only  the  zero  matrix 
but  any'  constant  times  the  identity.  Addition  of  any  of  these  "zero 
information"  matrices  leaves  the  corresponding  density  function, 
pQ,  unchanged  (after  normalizing). 

In  the  case  of  very  good  information  about  the  most  likely 
orientation,  or  small  error  covariances,  these  probability  density 
functions  ought  to  behave  like  the  multivariate  normal.  To  see  that 
this  is  the  case,  let  us  look  at  a  D  matrix  whose  density  function, 
PD,  is  concentrated  in  the  vicinity  of  b^,  the  group  identity. 

First,  bp  being  a  eigenvector  implies 

A 

0 

0 

0 


so,  by  symmetry  of  D, 


’  A  |  0  0  0  1 

|7T 

0  !  X 

L°  i 


(126) 


where  X  is  a  symmetric  3x3  matrix.  We  want  the  eigenvalues  of  X 
to  be  large  in  comparison  to  A;  as  with  all  such  D  matrices,  we  are 
free  at  this  point  to  subtract  AI  from  D,  the  effect  being  absorbed 
in  the  normalization  function,  k(D’),  where  D’  =  D-AI.  For  small 
variations  around  bn,  the  three  degrees  of  freedom  reside  in  the 
last  three  components  of  58,  as  we  have  seen  in  Eqs.  4  or  88.  Thus 


D’  =  D-AI  = 


"o  J  0  0  0 

~ol - 


(127) 
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does  indeed  look  like  a  multivariate  normal  density  in  the  three 
independent  differential  rotations  about  the  x,  y,  and  z  axes  under 
the  correspondence 


where  Y  =  E(00T),  6 


(128) 


‘WV’ 


and  E  is  the  expectation  operator. 


Noting  that  the  correspondence  in  Eq.  128  arises  from  the 

0  0 

approximation  sin  —  ~  y  for  small  0  in  the  definition  of  g  , 

0 

and  and  that  6q  involves  cos  y  ,  it  is  natural  to  ask  whether 

-IT  T 

D  might  be  obtainable  as  E(BB  ).  Certainly  the  matrix  E(BB  ) 

exists  and  is  symmetric  for  any  probability  density,  p^.  For  the 
uniform  density  function  Eq .  113,  it  is  not  difficult  to  evaluate 
E(66T), 


1 

2 


/ 

b|=i 


Pj(B)  88 


dg 


1/4 

1/4 

0 


0 

1/4 

1/4 


(129) 


where  dg  is  given  in  Eq.  105.  This  gives  D  =  41  and  shows  that  in 
this  case 


D_1  =  E(8BT) 


(130) 


for  an  appropriately  selected  D  from  the  class  (of  scalar  multiples 
of  I)  corresponding  to  the  uniform  density  function.  It  is  not 
difficult  to  show  that  each  class  J  D-l-X I ,  X  realj  ,  where  D  is 
symmetric,  has  one  choice  of  A  that  makes 

Trace  (D+AI)-1  =  1  .  (131) 
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The  trace  of  (D+AI)  is  the  sum  of  the  reciprocals  of  the  eigen¬ 
values  of  D+AI.  For  large  A  these  reciprocals  are  all  positive, 

but  small;  let  u .  be  the  eigenvalues  of  D  and  let  A  decrease  toward 
1 

the  largest  quantities  -tu,  and  the  trace  of  (D+AI)  will  increase, 
reaching  1  at  some  value  of  A. 

T 

The  significance  of  this  is  in  the  fact  that  E(8B  )  must 

have 


Trace  E(8BT)  =  1 


(132) 


because  68=1.  Is  it  the  case,  we  may  ask,  that  given  any  real 
symmetric  (4  x  4)  matrix  D 


E(66T) 


k(D) 


/bb’ 
I  8 1=1 


1  T 
j  6  DB 


dg  =  (D+AI) 


-1 


(133) 


Here  A  would  be  that  number  for  which  Eq .  131  is  satisfied.  No 
attempt  here  will  be  made  to  resolve  this  speculation. 

Consider  the  diagonalization  of  a  real  symmetric  (4  x  A) 
matrix  D  by  an  orthogonal  matrix  of  determinant  1: 


D  =  T(a,y)r  AT(a,v) 

where  T(a,B)  =  R(a)S(6)  as  in  Eq.  20,  and 

AQ  0  0  0 

0  A  0  0 

A  = 

0  0  A£  0 

_0  0  0  A3_ 


(134) 


(135) 


This  clearly  generalizes  Eqs.  119  and  120  to  two  types  of  "rota¬ 
tions"  that  arise  because  of  the  noncommutativity  of  composition 
of  rotations.  The  same  argument  used  to  establish  Eq.  121  shows 
that 
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The  density  function  Eq.  108  can  thus  be  reduced  by  a  change  of 

T 

variable  that  diagonalizes  the  quadratic  form,  8  D8,  to  the  function 


PA<«)  =  k<A)  e  4  (137) 

where  A  is  obtained  from  D  by  the  change  of  variables 

6  =  R (a)  S(y)B  .  (138) 

To  accomplish  this  reduction  to  standard  form  Eq.  137,  we  required 
not  one  but  two  transformations  of  the  stochastic  variable,  6,  in 
defining  the  new  stochastic  variable,  8.  By  commutativity  of  R  and 
S  matrices  (Eq.  16)  and  relationship,  Eq.  15,  this  transformation 
can  also  be  written 

S (6)  =  S (y)  S«3)  S(a)  ,  (139) 

where  the  nature  of  a  as  a  "precedent"  rotation  and  y  as  a  "subse¬ 
quent"  rotation  is  more  clearly  shown. 
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6.  ESTIMATION  OK  ORIENTATION 


The  fundamental  result  of  the  previous  section  that  we  will 
draw  on  is  the  means  of  calculating  the  most  likely  orientation, 

6*,  from  a  probability  density  function,  p^((3).  It  was  shown  that 

to  find  S*  from  D  it  was  necessary  to  look  for  the  eigenvector  cor¬ 
responding  to  the  smallest  eigenvalue.  This  answer  may  not  be 
unique;  two  eigenvalues  may  be  equal  or  nearly  so.  The  other  fea¬ 
ture  of  the  family  of  functions  chosen  is  that  the  product  of  two 
is  also  in  the  same  class  and  the  matrix  resulting  from  such  an  op¬ 
eration  is  the  sum  of  the  two  matrices  for  the  two  factors  (Eq.  123). 

So  far  we  have  seen  only  one  specific  example  of  a  D  matrix 
interpreted  physically:  D  =  41,  which  represents  "no  information" 
and  generates  the  uniform  probability  density.  Figures  4  and  5 
show  two  examples  of  the  shape  of  the  function 

1  T 

-  j  6  DB 

f(6)  =  e 

plotted  for  the  two-dimensional  analog  of  the  Euler  parameters. 

The  inner  curve  is  the  unit  circle,  | B |  =1,  and  the  outer  curve 
has  its  radius  increased  to  1  +  f(8).  In  Fig.  4,  the  eigenvectors 
of  D  are  aligned  with  the  coordinate  axes;  in  Fig.  5,  they  have 
been  rotated  as  indicated.  Following  the  discussion  of  the  matrix 
given  in  Eq .  126,  we  see  that  if 

10  0  0 

0  d2  0  0 

0  0  d2  0 

0  0  0  d3 

where  d^.d^.d^  denote  any  positive  real  values  very  much  larger 
than  1,  the  function  Pp(B)  is  a  maximum  at  6*  =  (1,0, 0,0).  A  ma¬ 
trix  D'  having  an  equally  well  defined  eigenvector  at  an  arbitrary 
point  3'  is  given  by 

D'  =  S(B')  DS(B')T  ,  (141) 
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T 

since  S(0')  carries  6'  into  (1,0*0, 0) .  There  is  another  choice  at 
hand  that  has  this  same  eigenvector: 

D"  =  R(0’)  DR(0')T  .  (142) 

The  distinction  between  D'  and  D"  when  0'  ^  (1000)  is  that  of 
the  order  in  which  0’  and  the  stochastic  variable  are  applied  as 
was  discussed  at  the  enc.  of  the  previous  section.  An  easy  way  to 
visualize  this  is  to  let  0'  become  a  function  of  time  with  0'(O)  = 
(1000).  Now 

S(0’)T  D'  S(0') 

equals  D  and  is  constant  in  time,  whereas 

D"*  (t)  =  S(0')T  D"S(0’)  =  E(0')T  D  1(0')  ,  (143) 

where  E(0')  is  defined  in  Eq.  17  and  D"’(t)  is  clearly  the  repre¬ 
sentation  of  the  constant  matrix  D  in  a  coordinate  frame  that  ro¬ 
tates  according  to  0 ' ( t ) ;  by  this  we  mean  that  the  choice  of  the 
three  coordinate  axes  in  Eqs.  1  and  4  is  a  rotating  frame,  as  is 
clear  from  Eq .  17.  Note  then  that  the  distinction  between  Eqs.  141 
and  142  must  disappear  in  the  special  case  where  d^  =  =  d^. 


Next,  consider  the  case  of  a  satellite  attitude  measuring 
system  such  as  a  vector  magnetometer  or  solar  attitude  deter  ..or 
that  provides  at  some  instant  of  time  the  information  that  a  body- 
fixed  unit  vector,  u,  is  aligned  with  an  inertial  unit  vector,  U. 
No  information  is  provided  about  the  orientation  of  the  spacecraft 
around  that  axis.  Given  only  these  data,  all  orientations  around 
this  axis  are  equally  likely.  Starting  with  any  a  such  that 


(144) 


(145) 


where  again  d  >>  1  (depending  on  the  precision  of  the  measuring 
instrument).  The  probability  density  function  pD,(0)  now  takes 

on  its  maximum  at  0  =  (1,0, 0,0),  but  this  point  is  not  unique. 
Since  any  0  of  the  form 
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is  carried  by  2(a)  into 


(146) 


which  is  also  an  eigenvector  of  the  diagonal  matrix  in  Eq.  145,  it 
follows  that  8*  is  an  eigenvector  of  D*  having  an  eigenvalue  of  2. 
Such  a  set  of  6'  describes  all  possible  orientations  obtained  by 
rotations  about  the  vector  u.  Since  the  vector  u  is  fixed  in  the 
body,  we  next  wish  to  precede  such  rotations  by  one  carrying  the 
inertial  vector  U  into  u;  call  this  vector  6. 

u  =  C(6)  U  .  (147) 


Writing 

D  =  R(6)  D’R(6)T  ,  (148) 

we  can  show  that  6  is  an  eigenvector  of  D.  Since  6  is  carried  by 
T 

R(6)  into  (1,0, 0,0),  which  we  have  shown  to  iv  an  eigenvector  of 
D'  with  eigenvalue  2  (the  smaller  of  2  and  d) ,  6  is  an  eigenvector 
of  D  and  Pp(8)  takes  on  an  extremum  (maximum)  at  6.  Each  of  the 

other  equally  likely  orientations  arises  from  a  S'  such  that 


8'  =  R(6)T6' 

(149) 

From  this  equation,  we  find 

S'  -  S(8')6 

(150) 

which  is  interpreted  as  S'  being  the  rotation  5  followed  by  the  ro¬ 
tation  8'  (Eq.  146). 
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Next  we  shall  give  a  numerical  example  to  show  how  two  such 
observations  may  be  combined  by  the  rule  of  adding  D-matrices  (mul¬ 
tiplying  density  functions)  to  yield  a  D-matrix  whose  eigenvector 
having  the  smallest  eigenvalue  provides  a  uinque  solution  of  the 
satellite  attitude.  Consider  first  the  simplest  case  in  which  the 
body-fixed  and  inertial  frames  happen  to  be  aligned;  the  Euler 
parameters  describing  this  relationship  are  (1,0, 0,0),  the  identity. 
Suppose  we  have  two  observations,  one  of  the  body  z-axis  and  one  of 
the  body  x-axis.  To  describe  the  first  by  a  D-matrix,  D^,  we  need 

to  find  a  and  6  according  to  Eqs.  144  and  147.  In  this  case  we  may 
choose  (neither  choice  is  unique)  a  =  6  =  (1,0, 0,0)  because  both  u 
and  U  are  the  respective  z-axis  unit  vectors,  (0,0,1),  of  their  co¬ 
ordinate  frames.  We  obtain 

2  0  0  0 

0  d  0  0 

0  0  d  0 

0  0  0  2 

and  will  reserve  the  choice  of  d. 

Next,  the  second  observation  is  defined  by  u  =  (1,0,0)  and 
U  =  (1,0,0),  and  we  assume  this  is  made  simultaneously  with  the 
first.  Here  again,  6  may  be  selected  as  (1,0, 0,0).  The  choice  of 
a  must  provide  the  point  transformation  (see  Appendix  A)  of  (1,0,0) 
Into  (0,0,1),  so  let 


(151) 


(152) 
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(154) 


Now  we  can  write  the  matrix  of  combined  information: 


D  -  n1  +  d2 


0  d+2 


0  2d 


(155) 


which,  for  large  d  has  a  unique  eigenvector,  (1,0, 0,0),  correspond¬ 
ing  to  its  smallest  eigenvalue  (4) . 

As  a  second  example,  we  let  the  orientation  of  the  space¬ 
craft  be  (all  calculations  carried  15  digits  although  we  show 
fewer  here) 


0.98255 

0.04971 

0.09942 

0.14913 


C(8*) 


/  0.93575  0.30293  -0.1804  \ 

-0.28316  0.95058  0.12733  1  . 

V 0.21019  -0.06803  0.97529/ 


We  pick  two  vectors  in  inertial  coordinates  as  observations 


-0.26726 
-0.53452 
.  0.80178 


0.96309 
-0.24077 
v  0.12039 
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These  are  mapped,  respectively,  into 

/-O. 55677  v 
ul  =  (  -0.33033  j 
V  0.76216  / 


/  0. 80654 \ 
u2  =  (-0.48626  J 
\  0.33622  / 

by  C(B*)  as  in  Eq.  147.  Now  treating  these  as  observations  with 


0,  we  reconstruct 
): 

8*  from  the 

sum  of  the  two 

D-matrices 

319.63 

-  712.53 

-1213.32 

-1046.39 

-  712.53 

8194.63 

-1929.06 

3249.87 

-1213.32 

-1929.06 

7917.01 

3360.56 

-1046.39 

3249.87 

3360.56 

3572.71 

t 

330.32 

-  112.05 

-1133.35 

\ 

-1370.29 

-  112.05 

1905.54 

3311.83 

-2104.12 

-1133.35 

3311.83 

8501.14 

697.31 

-1370.29 

-2104.12 

697.31 

9266.99 

D1+D2 


'  649.95 
824.59 
2346.67 
2416.69 


-  824.59 
10100.17 
1382.76 
1145.74 


-2346.67 

1382.76 

16418.16 

4057.88 


-2416.69' 

1154.74 

4057.88 

12839.70 


The  eigenvalues  and  eigenvectors  were  determined  by  the  Jacobi 
method  and  returned  6*  correct  to  at  least  ten  digits. 

When  observations  are  not  taken  simultaneously,  but  the 
motion  of  the  satellite  is  known  over  the  interval  between  observa¬ 
tions,  estimation  of  the  orientation  is  still  possible.  From 
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knowledge  of  w(t)  in  body-fixed  coordinates  we  can  produce  the 
transition  matrix 

•  (tj.tj)  =  S(f$2)  S(B1)T  (156) 

as  in  Eq.  31  without  actually  knowing  8^  or  &2,  which  represent  the 
spacecraft  orientation  at  t^  and  t2  when  the  two  observations  are 
taken.  The  estimate  of  6^  can  be  based  on 

D  =  $T  D2$  +  Dx  ,  (157) 

where  represents  the  observation  at  t^,  and  D2  that  at  t2> 

The  extension  of  this  method  to  the  case  in  which  u  is 
treated  as  a  stochastic  variable  is  not  addressed  here.  We  can 
speculate  that  the  convolution  of  two  density  functions  of  the  form 
Eq.  108  is  another  of  this  type,  and  further  that  if  w  assumes  a 
normal  density  that  a  differential  equation  for  propagation  of  D- 
matrices  would  represent  the  "diffusion"  of  these  probability  den¬ 
sities  and  a  sort  of  rotational  Brownian  motion  as  discussed  in 
Ref.  8  and  in  other  references  as  far  back  as  1928. 
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APPENDIX  A 

DEFINITION  OF  ROTATION  MATRICES 


The  confusion  that  steins  from  the  wide  variety  of  defini¬ 
tions  of  coordinate  frames,  Eulerian  angles,  and  rotation  matrices 
cannot  be  laid  aside  easily.  Reference  5,  for  example,  very  care¬ 
fully  discusses  the  freedom  of  choice  of  clockwise  versus  counter¬ 
clockwise  rotations  as  positive  rotations,  the  concept  of  "running 
coordinates,"  and  the  order  of  matrix  multiplication.  Consistent 
with  this,  Goldstein  writes  as  a  rotation  about  z  in  a  right-handed 
system  with  counterclockwise  rotations  being  positive  (see  Fig.  A-l) 

sin  <j>  0 

cos  <J>  0 

0  1 

He  completely  neglects  to  describe  this  as  an  "alias"  or  "coordinate 
transformation,"  and  the  casual  reader  might  be  surprised  that  for 
=  +90°  the  matrix  D  carries  (1,0,0)  into  (0,-l,0).  The  matrix 
C(B)  in  Eq.  6  is  similarly  defined  as  a  coordinate  transformation. 

Reference  3  provides  a  nice  discussion  of  "alias"  versus 
"alibi."  In  other  works,  notably  those  using  tensor  notation,  the 
terms  "covariant"  and  "contravariant"  appear.  The  choice  of  the 
"alias"  or  "coordinate  transformation"  or  "contravariant"  form 
above  is  necessary  if  the  matrix  products  are  to  be  taken  from 
right  to  left  in  "running  coordinates."  To  remember  this,  consider 
the  succession  of  two  rotations,  the  first  about  z  as  above,  and 
the  second  about  the  x-axis  after  it  has  been  repositioned  accord¬ 
ing  to  the  first  rotation.  We  choose  to  write  the  second  matrix 
in  "running  coordinates;"  that  is,  the  axis  of  x-rotation  is  to  be 
(1,0,0): 

j 1  0 

C  =»  10  cos  8 

\  0  -sin  9 


(cos  4> 
-sin  <J> 
0 


-  53  - 


THE  JOHNS  HOPKINS  UNIVERSITY 

APPLIED  PHYSICS  LABORATORY 

LAUREL  MARYLAND 


3,  3' 


Fig.  A-2  The  second  (subsequent)  rotation,  again  counter 
clockwise,  about  the  V  axis,  as  in  Ref.  5. 
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See  Fig.  A-2  for  a  description  of  this  relative  to  the  first  rota¬ 
tion.  Now  in  the  "alias"  interpretation  we  fix  a  point  whose  coor¬ 
dinates  are  initially  (cos  4>,  sin  0).  This  point  will  be  the 
final  x-axis  after  repositioning  the  coordinate  frame.  Indeed, 
applying  D  then  C  yields  (1  0  0). 


Using  the  "point  transformation"  or  "alibi"  method,  we  would 
piece  together  this  same  rotation  by  first  applying  a  rotation  to 
map 


and 


This  matrix 

is  DT. 

The 

/  cos 

♦\  / 

'  cos 

(  sin 

*  1  —  1 

sin 

The  subsequent  rotation  must  now  map 


and 


sin  0  sin  <f> 
-sin  9  cos  <j> 
cos  6 


and  is  clearly  not  CT.  It  is  in  fact  DTCTD,  so  the  product  is 


(DTCTD)  d1  -  dtct 


which  is  reasonably  enough  the  transpose  of  the  alias  representa¬ 
tion  CD.  Authors  who  prefer  the  alias  representation  probably  feel 
that  it  is  more  natural  because  the  matrices  in  running  coordinates 
(which  are  easiest  to  write)  are  applied  in  the  same  order  as  are 
the  successive  rotations.  Authors  who  prefer  the  alibi  probably 
feel  more  comfortable  watching  a  point  rotate  according  to  the  ri.ght- 
hand  rule,  e.g.,  (1 , 0 ,0)  — (cos  4> ,  sin  0)  for  "positive"  rotation 
about  z.  The  most  common  mistake  is  in  selecting  the  alibi  repre- 

T  T 

sentation  and  incorrectly  composing  two  rotations  as  C  D  (Ref.  6). 


THE  JOHNS  HOPKINS  UNIVERSITY 

APPLIED  PHYSICS  LABORATORY 

laurel  Maryland 


APPENDIX  B 


DEVELOPMENT  OK  e 


S(u>)l 


for  u>  CONSTANT 


“  0 

"“I 

~U2 

~u3  ~ 

“l 

0 

u3 

~u2 

“2 

"“3 

0 

U1 

f 

_  “3 

“2 

-U1 

0 

S(o>)  = 


by  definition,  which  it  is  convenient  to  partition,  letting 


(B.l) 


G  = 


3 

0 


-0), 


(B.2) 


into 


S(«o)  = 


0 

8 


_8 

G 


Murnaghan  (Ref.  7)  points  out  that  G  satisfies 
1 2 


G3  = 


G 

G 


(B.  3) 


(B.  4) 


and  that  therefore  e  is  expressible  as  a  quadratic  polynomial  in 
G.  It  is  natural  to  expect  the  same  of  S(w)  and  thus  of  S(to)t. 

We  observe  that 


Gg 


(B.5) 
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The  constant  term  is 


From  the  equality 


6(t)  =  S[e<t)]  bQ, 


and 

i  S(ai)t 
S[6(t)]  =  e 


(B. 11) 


(B.12) 
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